// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtModelFitFilter.C
// ************************************************************************* //

#include <avtModelFitFilter.h>
#ifdef PARALLEL
  #include<mpi.h>
  #include<BufferConnection.h>
#endif
#include<avtParallel.h>

#include<ParsingExprList.h>

int compare_floats (const void *, const void *);
int compare_doubles(const void *, const void *);

// ****************************************************************************
//  Method: avtModelFitFilter constructor
//
//  Programmer: rob -- generated by xml2avt
//  Creation:   Wed Feb 3 16:15:29 PST 2010
//
// ****************************************************************************

avtModelFitFilter::avtModelFitFilter()
{
    num_relats = 0;
}


// ****************************************************************************
//  Method: avtModelFitFilter destructor
//
//  Programmer: rob -- generated by xml2avt
//  Creation:   Wed Feb 3 16:15:29 PST 2010
//
//  Modifications:
//
// ****************************************************************************

avtModelFitFilter::~avtModelFitFilter()
{
}

// ****************************************************************************
//  Method:  avtModelFitFilter::Create
//
//  Programmer: rob -- generated by xml2avt
//  Creation:   Wed Feb 3 16:15:29 PST 2010
//
// ****************************************************************************

avtFilter *
avtModelFitFilter::Create()
{
    return new avtModelFitFilter();
}


// ****************************************************************************
//  Method:      avtModelFitFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: rob -- generated by xml2avt
//  Creation:   Wed Feb 3 16:15:29 PST 2010
//
// ****************************************************************************

void
avtModelFitFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const ModelFitAtts*)a;
}


// ****************************************************************************
//  Method: avtModelFitFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtModelFitFilter with the given
//      parameters would result in an equivalent avtModelFitFilter.
//
//  Programmer: rob -- generated by xml2avt
//  Creation:   Wed Feb 3 16:15:29 PST 2010
//
// ****************************************************************************

bool
avtModelFitFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(ModelFitAtts*)a);
}


// ****************************************************************************
//  Method: avtModelFitFilter::xecuteData
//
//  Purpose:
//      Sends the specified input and output through the ModelFit filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: rob -- generated by xml2avt
//  Creation:   Wed Feb 3 16:15:29 PST 2010
//
//  Modifications:
//    Eric Brugger, Wed Jul 30 19:34:44 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************
#include<vtkDataArray.h>
#include<vtkDoubleArray.h>
#include<vtkFloatArray.h>
#include<vtkRectilinearGrid.h>

#include<avtCallback.h>
#include<avtDatabase.h>
#include<avtDatabaseMetaData.h>

avtDataRepresentation *
avtModelFitFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    intVector numTups = atts.GetNumTups();
    intVector numVars = atts.GetNumVars();
    stringVector Vars = atts.GetVars();

    if(!num_relats)
    for(size_t i = 0; i < numTups.size(); i++)
        if(numTups[i] > 0)
        num_relats++;
    
    int tup_index, var_index;
    vtkDataArray *darray;
    char secVarName[1024];
     
    //outer loop grabs first correctly specified variable
    //innards for determining centering
    for(int i = 0; i < num_relats; i++)
    {
        if(numTups[i] < 1)
            continue;

        tup_index = 0;
        var_index = 0;
        for(int j = 0; j < i; j++)
        {
            tup_index += numTups[j]*numVars[j];
            var_index += numVars[j];
        }
    
        sprintf(secVarName, "%s", Vars[var_index].c_str());                

        darray = in_ds->GetPointData()->GetScalars(secVarName);                            
        if(darray == NULL)
            Centering = AVT_ZONECENT;
        else
            Centering = AVT_NODECENT;
        break;
    }
    
    vtkFloatArray *arr;
    int numValues;
    std::string varname;

    //Creating new variable for storing model number
    if(in_ds->GetPointData()->GetScalars("operators/ModelFit/model") == NULL &&
        in_ds->GetCellData()->GetScalars("operators/ModelFit/model") == NULL)
    {
        arr = vtkFloatArray::New();

        if(Centering == AVT_ZONECENT)
            numValues = in_ds->GetNumberOfCells();
        else 
            numValues = in_ds->GetNumberOfPoints();
    
        arr->SetNumberOfComponents(1);
        arr->SetNumberOfValues(numValues);
    
        varname = "operators/ModelFit/model"; 
        arr->SetName(varname.c_str());
    
        if(Centering == AVT_ZONECENT)
        {
            in_ds->GetCellData()->AddArray(arr);
            in_ds->GetCellData()->SetActiveScalars(varname.c_str());
        }
        else
        {
            in_ds->GetPointData()->AddArray(arr);
            in_ds->GetPointData()->SetActiveScalars(varname.c_str());
        }
    }

    vtkFloatArray *arr2;

    //Creating variable for storing distance
    if(in_ds->GetPointData()->GetScalars("operators/ModelFit/distance") == NULL &&
        in_ds->GetCellData()->GetScalars("operators/ModelFit/distance") == NULL)
    {
        arr2 = vtkFloatArray::New();
    
        if(Centering == AVT_ZONECENT)
            numValues = in_ds->GetNumberOfCells();
        else 
            numValues = in_ds->GetNumberOfPoints();
    
        arr2->SetNumberOfComponents(1);
        arr2->SetNumberOfValues(numValues);

        varname = "operators/ModelFit/distance"; 
        arr2->SetName(varname.c_str());
    
        if(Centering == AVT_ZONECENT)
        {
            in_ds->GetCellData()->AddArray(arr2);
            in_ds->GetCellData()->SetActiveScalars(varname.c_str());
        }
        else
        {
            in_ds->GetPointData()->AddArray(arr2);
            in_ds->GetPointData()->SetActiveScalars(varname.c_str());
        }
    }

    out_ds.push_back(in_ds);
    return in_dr;
}

//    Mark C. Miller, Wed Aug 22 19:57:54 PDT 2012
//    Removed superfolous code allocating and initializing local variables,
//    score_mins and score_maxs. Nothing in the function was using these
//    variables.  
void
avtModelFitFilter::PostExecute()
{  
    intVector numTups = atts.GetNumTups();
    intVector numVars = atts.GetNumVars();
    stringVector Vars = atts.GetVars();
    doubleVector Tuples = atts.GetTuples();
    doubleVector thold  = atts.GetThold();
    unsignedCharVector StatTuples = atts.GetStatTuples();
    intVector tup_match;
    doubleVector mins;
    doubleVector point;

    int tup_index, var_index;
    vtkDataArray *darray;
    char secVarName[1024];
    int numValues;
    size_t varFinder;

    double distance;
    float *scalars;
    bool good_match;

    int relat_index;

    double opt_distance;

    intVector selectionType = atts.GetSelectionType();
    intVector distanceType  = atts.GetDistanceType();
    intVector inputSpace    = atts.GetInputSpace(); 


    //Since there may be overlap between models in variables specified,
    //create a set of all unique variables
    for(size_t i = 0; i < Vars.size(); i++)
    {
        sprintf(secVarName, "%s", Vars[i].c_str());
        createVS(secVarName);
    }
    
    //determine for each variable whether or not statistics or 
    //histograms need to be calculated
    for(int i = 0, currentVar = 0, num_to_remove = 0; i < num_relats; i++)
    {
        if(selectionType[i] == 1 || inputSpace[i] == 1)
            for(int j = 0; j < numVars[i]; j++, currentVar++, num_to_remove++)
            {
                for(varFinder = 0; varFinder < VS.size(); varFinder++)
                    if(!strcmp(VS[varFinder]->name, Vars[currentVar].c_str()))
                        break;
                VS[varFinder]->calculate_stats = true;
            }
        currentVar -= num_to_remove;
        num_to_remove = 0;

        if(selectionType[i] == 3 || inputSpace[i] == 3)
            for(int j = 0; j < numVars[i]; j++, currentVar++, num_to_remove++)
            {
                for(varFinder = 0; varFinder < VS.size(); varFinder++)
                    if(!strcmp(VS[varFinder]->name, Vars[currentVar].c_str()))
                        break;
                VS[varFinder]->calculate_hists = true;
            }
        currentVar -= num_to_remove;
        num_to_remove = 0;

        currentVar += numVars[i];
    }
    
    //also force stats calculation if stat shortcuts specified
    for(int i = 0; i < num_relats; i++)
    {
        tup_index = 0;
        var_index = 0;
        for(int j = 0; j < i; j++)
        {
            tup_index += numTups[j]*numVars[j];
            var_index += numVars[j];
        }
        for(int j = 0; j < numTups[i]*numVars[i]; j++)
            if(StatTuples[tup_index+j] == 'a' || StatTuples[tup_index+j] == 'M' || 
                    StatTuples[tup_index+j] == 'm')
            {
                for(varFinder = 0; varFinder < VS.size(); varFinder++)
                    if(!strcmp(VS[varFinder]->name, Vars[var_index + (j%numVars[i])].c_str()))
                        break;

                VS[varFinder]->calculate_stats = true;
            }
    }
    
    //Go through all variables specified and calculate
    //statistics for each.
    calculateVariableStats();
    
#ifdef PARALLEL
    int numProcs = 1;
    int myRank   = 0;
    
    MPI_Comm_rank(VISIT_MPI_COMM, &myRank);
    MPI_Comm_size(VISIT_MPI_COMM, &numProcs);

    double     *sum;
    double     *max;
    double     *min;
    double     *non_zero_min;
    int *num_points;
    int totVars  = VS.size();
    int number_on = 0;

    double *temp_hist;
    double *merged_ranges;
    
    num_points = (int    *)malloc(sizeof(int)   *totVars);
    sum        = (double *)malloc(sizeof(double)*totVars);
    max        = (double *)malloc(sizeof(double)*totVars);
    min        = (double *)malloc(sizeof(double)*totVars);
    non_zero_min = (double *)malloc(sizeof(double)*totVars);

    merged_ranges = (double *)malloc(sizeof(double)*(BINS*numProcs));
    temp_hist = (double *)malloc(sizeof(double)*BINS);

    for(int i = 0; i < totVars; i++)
    {
        if(VS[i]->calculate_stats)
        {
            num_points[i] = VS[i]->num_points;
            sum[i]        = VS[i]->sum;
            max[i]        = VS[i]->max;
            min[i]        = VS[i]->min;
            non_zero_min[i] = VS[i]->non_zero_min; 
    
            MPI_Allreduce(&num_points[i], &(VS[i]->num_points), 1, MPI_INT, MPI_SUM, VISIT_MPI_COMM);
            MPI_Allreduce(&sum[i], &(VS[i]->sum), 1, MPI_DOUBLE, MPI_SUM, VISIT_MPI_COMM);
            MPI_Allreduce(&max[i], &(VS[i]->max), 1, MPI_DOUBLE, MPI_MAX, VISIT_MPI_COMM);
            MPI_Allreduce(&min[i], &(VS[i]->min), 1, MPI_DOUBLE, MPI_MIN, VISIT_MPI_COMM);
            MPI_Allreduce(&non_zero_min[i], &(VS[i]->non_zero_min), 1, MPI_DOUBLE, MPI_MIN, VISIT_MPI_COMM);
        }

        if(VS[i]->calculate_hists)
        {
            for(int np = 0; np < numProcs; np++)
            {
                if(np == myRank)
                    memcpy(temp_hist, VS[i]->hist.range_ends, BINS*sizeof(double));
            
                MPI_Bcast(temp_hist, BINS, MPI_DOUBLE, i, VISIT_MPI_COMM);
           
                for(int j = 0; j < BINS; j++)
                    merged_ranges[number_on++] = temp_hist[j];
            }
            
            VS[i]->hist.count      = VS[i]->num_points / BINS;
            VS[i]->hist.mod_number = VS[i]->num_points % BINS;
        
            qsort(merged_ranges, BINS*numProcs, sizeof(double), compare_doubles);
        
            for(int j = 0, k = 0; j < numProcs*BINS; j += numProcs, k++)
                VS[i]->hist.range_ends[k] = merged_ranges[j];
            number_on = 0;
        }
    }

    delete num_points;
    delete sum;
    delete max;
    delete min;
    delete non_zero_min;
    delete(temp_hist);
    delete(merged_ranges);
    
#endif

    for(size_t i = 0; i < VS.size(); i++)
        if(VS[i]->num_points)
            VS[i]->average = VS[i]->sum/VS[i]->num_points;

    /*
    puts("-------before------");
    for(int i = 0; i < Tuples.size(); i++)
        printf("Tuples[%d]: %lf\n", i, Tuples[i]);
        */

    //do space converts on thresholds and tuples
    //for tuples, stick stats in if needed
    for(int i = 0; i < num_relats; i++)
    {
        tup_index = 0;
        var_index = 0;
        for(int j = 0; j < i; j++)
        {
            tup_index += numTups[j]*numVars[j];
            var_index += numVars[j];
        }
        
        for(int j = 0; j < numTups[i]*numVars[i]; j++)
        {
            sprintf(secVarName, "%s", Vars[var_index+(j%numVars[i])].c_str());
            for(varFinder = 0; varFinder < VS.size(); varFinder++)
                if(!strcmp(VS[varFinder]->name, secVarName))
                    break;
 
            //we have a thold specified
            if(thold[tup_index+j] > 0)
                if(inputSpace[i] != selectionType[i])
                {
                    thold[tup_index+j] = convertToVariableSpace(thold[tup_index+j], varFinder, inputSpace[i]);
                    thold[tup_index+j] = spaceConvert(thold[tup_index+j], varFinder, selectionType[i]);
                }

            //replace Tuples with stats
            if(StatTuples[tup_index+j] == 'a' || StatTuples[tup_index+j] == 'M' || 
                    StatTuples[tup_index+j] == 'm')
            {
                if(StatTuples[tup_index+j] == 'a')
                    Tuples[tup_index+j] = VS[varFinder]->average;
                else if (StatTuples[tup_index+j] == 'M')
                    Tuples[tup_index+j] = VS[varFinder]->max;
                else if (StatTuples[tup_index+j] == 'm')
                    Tuples[tup_index+j] = VS[varFinder]->min;

                //printf("Tuples[%d] before: %lf\n", tup_index+j, Tuples[tup_index+j]);
                Tuples[tup_index+j] = spaceConvert(Tuples[tup_index+j], varFinder, selectionType[i]);
                //printf("Tuples[%d] after: %lf\n", tup_index+j, Tuples[tup_index+j]);
            }
            else if(inputSpace[i] != selectionType[i])
            {
                //printf("Tuples[%d] before: %lf\n", tup_index+j, Tuples[tup_index+j]);
                Tuples[tup_index+j] = convertToVariableSpace(Tuples[tup_index+j], varFinder, inputSpace[i]);
                //printf("Tuples[%d] after cTVS: %lf\n", tup_index+j, Tuples[tup_index+j]);
                Tuples[tup_index+j] = spaceConvert(Tuples[tup_index+j], varFinder, selectionType[i]);
                //printf("Tuples[%d] after SC: %lf\n", tup_index+j, Tuples[tup_index+j]);
                //printf("Used %lf for min, %lf for max\n", VS[varFinder]->min, VS[varFinder]->max);
            }
        }
    }
    /*
    puts("-------after------");
    for(int i = 0; i < Tuples.size(); i++)
        printf("Tuples[%d]: %lf\n", i, Tuples[i]);
        */
    atts.SetTuples(Tuples);
    atts.SetThold(thold);

    //
    for(size_t ds = 0; ds < out_ds.size(); ds++){
        numValues = 0;
       
        if(Centering == AVT_NODECENT)
            numValues = out_ds[ds]->GetNumberOfPoints();
        else
            numValues = out_ds[ds]->GetNumberOfCells();
    
        for(int p = 0; p < numValues; p++)
        {
            opt_distance = -1; //added for
            relat_index  = -1; //distance use
        
            for(int i = 0; i < num_relats; i++)
            { //added for distances
                if(numTups[i] < 1)
                    continue;
        
                tup_index = 0;
                var_index = 0;
        
                for(int j = 0; j < i; j++)
                {
                    tup_index += numTups[j]*numVars[j];
                    var_index += numVars[j];
                }
        
                grabOnePoint(&point, numVars[i], var_index, selectionType[i], p, out_ds[ds]->GetPointData(), 
                        out_ds[ds]->GetCellData());

                findMatches(&mins, &tup_match, point, numVars[i], numTups[i], tup_index, var_index);
        
                //Find winning relationship
                if(tup_match.size() > 0)
                {
                    distance = calculateDistance(mins, distanceType[i]);
    
                    good_match = true;
                    for(size_t j = 0; j < mins.size(); j++)
                        if(mins[j] < 0)
                            good_match = false;
            
                    if(good_match)
                        if(distance < opt_distance || opt_distance < 0)
                        {
                            opt_distance = distance;
                            relat_index = i;
                        }
                }
            }
        
            if(num_relats && out_ds.size() && p < numValues){
                if(Centering == AVT_NODECENT)
                    darray = out_ds[ds]->GetPointData()->GetScalars("operators/ModelFit/model");
                else
                    darray = out_ds[ds]->GetCellData()->GetScalars("operators/ModelFit/model");
        
                scalars = (float *)darray->GetVoidPointer(0);
                scalars[p] = relat_index+1;
    
                if(Centering == AVT_NODECENT)
                    darray = out_ds[ds]->GetPointData()->GetScalars("operators/ModelFit/distance");
                else
                    darray = out_ds[ds]->GetCellData()->GetScalars("operators/ModelFit/distance");
        
                scalars = (float *)darray->GetVoidPointer(0);
                if(opt_distance > 0)
                    scalars[p] = opt_distance;
                else
                    scalars[p] = -1;
            }
        }
    }
}


avtContract_p
avtModelFitFilter::ModifyContract(avtContract_p in_spec){
    std::string db = GetInput()->GetInfo().GetAttributes().GetFullDBName();
    ref_ptr<avtDatabase> dbp = avtCallback::GetDatabase(db, 0, NULL);
    avtDatabaseMetaData *md = dbp->GetMetaData(0);

    numTimesteps = (int)md->GetTimes().size();
    activeTs = in_spec->GetDataRequest()->GetTimestep();

    std::string meshName = GetInput()->GetInfo().GetAttributes().GetMeshname();

    pipelineVar = in_spec->GetDataRequest()->GetVariable();
    new_pipelineVar = strdup(pipelineVar);

    const stringVector curListedVars = atts.GetVars();

    if(!curListedVars.size())
    return in_spec;
  
    if((!strcmp(pipelineVar, "operators/ModelFit/model")) || 
       (!strcmp(pipelineVar, "operators/ModelFit/distance")))
        for(int i = 0; ; i++)
          if((strcmp(curListedVars[i].c_str(), "operators/ModelFit/model")) && 
             (strcmp(curListedVars[i].c_str(), "operators/ModelFit/distance")))
            {
                strcpy(new_pipelineVar, curListedVars[i].c_str());
                break;
            }

    avtDataRequest_p aDR  = new avtDataRequest(in_spec->GetDataRequest(), new_pipelineVar);
    aDR->SetOriginalVariable(pipelineVar);
    avtContract_p outSpec = new avtContract(in_spec, aDR);

    const char *curListedVar;
    std::vector<CharStrRef> curSecondaryVars = outSpec->GetDataRequest()->GetSecondaryVariables();
    size_t listedVarNum, secVarNum;
    char secVarName[1024];

    for(listedVarNum = 0; listedVarNum < curListedVars.size(); listedVarNum++){
        curListedVar = curListedVars[listedVarNum].c_str();
    
        if(strcmp(curListedVar, pipelineVar)){
            for(secVarNum = 0; secVarNum < curSecondaryVars.size(); secVarNum++)
                if(!strcmp(*curSecondaryVars[secVarNum], curListedVar))
                    break;
            if(secVarNum >= curSecondaryVars.size())
                outSpec->GetDataRequest()->AddSecondaryVariable(curListedVar);
            sprintf(secVarName, "%s", curListedVar);
        }
        curSecondaryVars = outSpec->GetDataRequest()->GetSecondaryVariables();
    }
    
    curSecondaryVars = outSpec->GetDataRequest()->GetSecondaryVariables();
    for(secVarNum = 0; secVarNum < curSecondaryVars.size(); secVarNum++){
        if(!strcmp(*curSecondaryVars[secVarNum], "operators/ModelFit/distance"))
            outSpec->GetDataRequest()->RemoveSecondaryVariable("operators/ModelFit/distance");

        if(!strcmp(*curSecondaryVars[secVarNum], "operators/ModelFit/model"))
            outSpec->GetDataRequest()->RemoveSecondaryVariable("operators/ModelFit/model");
    }
    return outSpec;
}

void
avtModelFitFilter::UpdateDataObjectInfo()
{
    avtDataAttributes &outatts = GetOutput()->GetInfo().GetAttributes();

    outatts.AddVariable("operators/ModelFit/model");
    outatts.ValidVariable("operators/ModelFit/model");
    outatts.SetVariableDimension(1, "operators/ModelFit/model");
    outatts.SetCentering(AVT_NODECENT, "operators/ModelFit/model");
    outatts.SetVariableType(AVT_SCALAR_VAR, "operators/ModelFit/model");

    outatts.AddVariable("operators/ModelFit/distance");
    outatts.ValidVariable("operators/ModelFit/distance");
    outatts.SetVariableDimension(1, "operators/ModelFit/distance");
    outatts.SetCentering(AVT_NODECENT, "operators/ModelFit/distance");
    outatts.SetVariableType(AVT_SCALAR_VAR, "operators/ModelFit/distance");
}

int compare_floats(const void *a, const void *b)
{
    float *arg1 = (float *)a;
    float *arg2 = (float *)b;
    if(*arg1 == *arg2) return 0;
    if(*arg1 < *arg2) return -1;
    return 1;
}

int compare_doubles(const void *a, const void *b)
{
    double *arg1 = (double *)a;
    double *arg2 = (double *)b;
    if(*arg1 == *arg2) return 0;
    if(*arg1 < *arg2) return -1;
    return 1;
}

#include<limits>
void
avtModelFitFilter::createVS(char *secVarName)
{
    size_t varFinder;
    variableStats *vs;

    for(varFinder = 0; varFinder < VS.size(); varFinder++)
        if(!strcmp(VS[varFinder]->name, secVarName))
            break;
    if(varFinder != VS.size())
        return;

    vs = (variableStats *)malloc(sizeof(variableStats));
    vs->name = strdup(secVarName);
    vs->sum = 0;
    vs->max = 0;
    vs->min = std::numeric_limits<double>::max();
    vs->non_zero_min = std::numeric_limits<double>::max();
    vs->num_points = 0;
    vs->average = 0;
    vs->calculate_stats = false;
    vs->calculate_hists = false;
    vs->hist.range_ends = (double *)malloc(sizeof(double)*BINS);
    VS.push_back(vs);
}

void
avtModelFitFilter::calculateVariableStats()
{
    int numValues;
    int single_count;
    int single_mod_number;
    int number_on;
    int total_num_values;

    float  *scalars;
    float  *sorters;
    double *merged_ranges;

    vtkDataArray *darray;

    //loop over VS
    for(size_t varFinder = 0; varFinder < VS.size(); varFinder++)
    {
        merged_ranges = (double *)malloc(sizeof(double)*(out_ds.size()*BINS));
        total_num_values = 0;
        number_on = 0;

        for(size_t ds = 0; ds < out_ds.size(); ds++)
        {
            if(Centering == AVT_NODECENT)                                            
            {
                numValues = out_ds[ds]->GetNumberOfPoints();
                darray = out_ds[ds]->GetPointData()->GetScalars(VS[varFinder]->name);
            }
            else                                                                     
            {
                numValues = out_ds[ds]->GetNumberOfCells();
                darray = out_ds[ds]->GetCellData()->GetScalars(VS[varFinder]->name); 
            }
        
            if(darray == NULL)                                                       
            {                   
                // TODO: We need to throw an exception or something sensible here
                printf("Variable errors (possibly centering) -- bailing\n");         
                //exit(1);                                 
            }                                        
        
            total_num_values += numValues;
            scalars = (float *)darray->GetVoidPointer(0);                             
        
            if(VS[varFinder]->calculate_stats)
            {
                for(int i = 0; i < numValues; i++)
                {
                    if(!VS[varFinder]->num_points)
                    {
                        VS[varFinder]->min = scalars[i];
                        VS[varFinder]->max = scalars[i];
                        VS[varFinder]->num_points = 1;
                    }
                    else
                    {
                        if(scalars[i] > VS[varFinder]->max)
                            VS[varFinder]->max = scalars[i];
                        else if(scalars[i] < VS[varFinder]->min)
                            VS[varFinder]->min = scalars[i];
                        if(scalars[i] < VS[varFinder]->non_zero_min && scalars[i])
                            VS[varFinder]->non_zero_min = scalars[i];
                    }
                    VS[varFinder]->sum+=scalars[i];
                }
            }
        
            if(VS[varFinder]->calculate_hists)
            {
                sorters = (float *)malloc(sizeof(float)*numValues);
                memcpy(sorters, scalars, numValues*sizeof(float));
    
                qsort(sorters, numValues, sizeof(float), compare_floats);
    
                single_count = numValues/BINS;
                single_mod_number = numValues % BINS;
                for(int i = single_count-1, j = 0; i < numValues; i+= single_count, j++)
                {
                    if(j < single_mod_number)
                        i++;
                    merged_ranges[number_on] = sorters[i];
                    number_on++;
                }
                delete(sorters);
            }
            VS[varFinder]->num_points = total_num_values;
        }
    
        
        if(VS[varFinder]->calculate_hists)
        {
            VS[varFinder]->hist.count      = total_num_values / BINS;
            VS[varFinder]->hist.mod_number = total_num_values % BINS;
        
            qsort(merged_ranges, out_ds.size()*BINS, sizeof(double), compare_doubles);
        
            for(size_t i = 0, j = 0; i < out_ds.size()*BINS; i+= out_ds.size(), j++)
                VS[varFinder]->hist.range_ends[j] = merged_ranges[i];
        }
        delete merged_ranges;
    }
}

void 
avtModelFitFilter::grabOnePoint(doubleVector *point, int numVars, int var_index, int selectionType, int which_point, 
    vtkPointData *inPointData, vtkCellData *inCellData) 
{
    stringVector Vars = atts.GetVars();
    float *scalars;
    vtkDataArray *darray;
    char secVarName[1024];

    for(int i = (int)point->size()-1; i >= 0; i--)
        point->erase(point->begin()+i); 

    size_t varFinder;

    for(int j = 0; j < numVars; j++)
    {//grabbing vars in a relationship - 1111111111111111111111111111111111111111
        sprintf(secVarName, "%s", Vars[var_index+j].c_str());                //11
                                                                             //11
        if(Centering == AVT_NODECENT)                                        //11
            darray = inPointData->GetScalars(secVarName);                    //11
        else                                                                 //11
            darray = inCellData->GetScalars(secVarName);                     //11
                                                                             //11
        if(darray == NULL){                                                  //11
            printf("Variable errors (possibly centering) -- bailing\n");     //11
            //exit(1);                                                         //11
        }                                                                    //11
        scalars= (float *)darray->GetVoidPointer(0);                         //11
                                                                             //11
        for(varFinder = 0; varFinder < VS.size(); varFinder++)               //11
        if(!strcmp(VS[varFinder]->name, secVarName))                         //11
            break;                                                           //11
        point->push_back(spaceConvert(scalars[which_point],                  //11
                    varFinder, selectionType));                              //11
    }//11111111111111111111111111111111111111111111111111111111111111111111111111

    for(int i = (int)Vars.size()-1; i >=0; i--)
        Vars.erase(Vars.begin()+i);
}

void
avtModelFitFilter::findMatches(doubleVector *mins, intVector *tup_match1, doubleVector point, 
    int numVars, int numTups, int tup_index, int var_index)
{
    doubleVector Tuples = atts.GetTuples();
    doubleVector thold  = atts.GetThold();
    stringVector Vars = atts.GetVars();
    intVector tup_match2;
    double min;
    int l;
    size_t varFinder;
    double delta;

    // Is the relationship a fit at all? /////////////////
    for(int i = (int)mins->size()-1; i >= 0; i--)
        mins->erase(mins->begin()+i); 

    for(int j = 0; j < numVars; j++)
    {
        for(varFinder = 0; varFinder < VS.size(); varFinder++)
            if(!strcmp(VS[varFinder]->name, Vars[var_index+j].c_str()))
                break;
    
        //Here we get from which tuple has a value closest to a point
        min = fabs(point[j] - Tuples[tup_index+j]);
    
        for(int k = 1; k < numTups; k++)
        {
            delta = fabs(point[j] - Tuples[tup_index+k*numVars+j]);
            if(delta < min)
                min = delta;
        }
    
        mins->push_back(min);
    }

    //find all tuples which have a min for the first var in a relationship
    for(int i = (int)tup_match1->size()-1; i >= 0; i--)
        tup_match1->erase(tup_match1->begin()+i);

    for(int k = 0; k < numTups; k++)
    {
        for(varFinder = 0; varFinder < VS.size(); varFinder++)
        if(!strcmp(VS[varFinder]->name, Vars[var_index+0].c_str()))
            break;

        delta = fabs(point[0] - Tuples[tup_index+k*numVars+0]);
        if(delta <= thold[tup_index+k*numVars+0] || thold[tup_index+k*numVars+0] <= 0)
            if(delta == (*mins)[0])
                tup_match1->push_back(k);
    }

    //now, find all tuples which match the rest of the variables
    for(int j = 1; j < numVars; j++)
    {
        for(varFinder = 0; varFinder < VS.size(); varFinder++)
            if(!strcmp(VS[varFinder]->name, Vars[var_index+j].c_str()))
                break;

        for(int k = 0; k < numTups; k++)
        {
            delta = fabs(point[j] - Tuples[tup_index+k*numVars+j]);
            if(delta <= thold[tup_index+k*numVars+0] || thold[tup_index+k*numVars+0] <= 0)
                if(delta == (*mins)[j])
                    tup_match2.push_back(k);
        }
    
        //make sure the min value comes from the same tuple for each variable
        for(int k = (int)tup_match1->size() - 1; k >= 0; k--)
        {
            for(l = 0; l < (int)tup_match2.size(); l++)
                if((*tup_match1)[k] == tup_match2[l])
                    break;
            if(l == (int)tup_match2.size())
                tup_match1->erase(tup_match1->begin()+k);
        }
    }
    for(int i = (int)tup_match2.size()-1; i >= 0; i--)
        tup_match2.erase(tup_match2.begin()+i); 
}


double 
avtModelFitFilter::calculateDistance(doubleVector mins, int distanceType)
{
    double distance = 0;
    if(!distanceType)
    {/* Euclidean Distance ////////////////////////ED*/
        for(size_t j = 0; j < mins.size(); j++)      //ED
        distance += pow(mins[j], 2);              //ED
        return sqrt(distance);                    //ED*/
    }
    if(distanceType == 1) 
    {/* Manhattan Distance /////////////////////////MD*/
    for(size_t j = 0; j < mins.size(); j++)          //MD
        distance += fabs(mins[j]);                //MD
    return distance;                              //MD*/
    }
    if(distanceType == 2)
    {/* Chebyshev Distance /////////////////////////CD*/
    double max = mins[0];                         //CD
    for(size_t j = 1; j < mins.size(); j++)          //CD
        if(mins[j] > max)                         //CD
        max = mins[j];                            //CD
    return max;                                   //CD*/
    }
    return distance; /// TODO: check fix for return of non-void warning
}

double avtModelFitFilter::convertToVariableSpace(double value, size_t varFinder, int fromType)
{
    if(!fromType)
        return value;

    if(fromType == 1)  //0-1
        return (value * (VS[varFinder]->max - VS[varFinder]->min) + 
            VS[varFinder]->min);
            
    if(fromType == 2)
        return (pow(2, value));

    int i;
    double number_below = value * VS[varFinder]->num_points;

    for(i = 0; i < BINS; i++)
    {
        if(number_below > 0)
        {
            number_below -= VS[varFinder]->hist.count;
            if(i < VS[varFinder]->hist.mod_number)
                number_below -= 1;
        }
        else
            break;
    }
    i--;

    number_below += VS[varFinder]->hist.count;
    if(i < VS[varFinder]->hist.mod_number)
        number_below += 1;

    if(i < VS[varFinder]->hist.mod_number)
        number_below /= (VS[varFinder]->hist.count+1);
    else
        number_below /= VS[varFinder]->hist.count;

    if(i < 2)
        return number_below * VS[varFinder]->hist.range_ends[0];
    return VS[varFinder]->hist.range_ends[i-1] + 
        (number_below * (VS[varFinder]->hist.range_ends[i] - VS[varFinder]->hist.range_ends[i-1]));

}
double 
avtModelFitFilter::spaceConvert(double value, size_t varFinder,  int selectType)
{
    if(!selectType)
        return value;

    if(selectType == 1)  //0-1)
    return (value-VS[varFinder]->min)/
        (VS[varFinder]->max-VS[varFinder]->min);

    if(selectType == 2)  //log
    {
        const double log2inv = 1.0 / log(2.0);
        if(value)
            return log(value) * log2inv;
        return log(VS[varFinder]->non_zero_min) * log2inv;
    }
    
    //probablity
    int number_below = 0;
    int i;
    double percent_o_bin;

    for(i = 0; i < BINS; i++)
    {
        if(VS[varFinder]->hist.range_ends[i] <= value)
        {
            number_below += VS[varFinder]->hist.count;
            if(i < VS[varFinder]->hist.mod_number)
                number_below++;
        }
        else
            break;
    }

    if(i == BINS)
        return 1;

    percent_o_bin = (value-VS[varFinder]->hist.range_ends[i])/
        (VS[varFinder]->hist.range_ends[i+1]-VS[varFinder]->hist.range_ends[i]);
    
    if(percent_o_bin < 0)
        percent_o_bin = 0;
    if(i < VS[varFinder]->hist.mod_number)
        number_below += (int)(percent_o_bin * (VS[varFinder]->hist.count+1));
    else
    
        number_below += (int)(percent_o_bin * VS[varFinder]->hist.count);

    return (double)((number_below*1.0)/(VS[varFinder]->num_points*1.0));
}
